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Abstract. 

The phase transitions at finite temperatures in the systems described by the 
Bose-Fermi-Hubbard model are investigated in this work in the framework of the 
selfconsistent random phase approximation. The case of the hard-core bosons is 
considered and the pseudospin formalism is used. The density-density correlator is 
calculated in the random phase approximation and the possibilities of transitions 
from superfluid to supcrsolid phases are investigated. It is shown that the transitions 
between uniform and charge ordered phases can be of the second or the first order, 
depending on the system parameters. 
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1. Introduction 

Properties of the systems of ultracold atomic gases confined in optical lattices are 
intensively studied in the last years both theoretically and experimentally [H [2j [3J 111 E] . 
Special attention is paid to the mixture of bosons and spin polarised fermions (e.g., 6 Li- 
7 Li, 40 K- 87 Rb, 6 Li- 87 Rb atoms). Such systems can be well described by the Bose-Fermi- 
Hubbard model (BFHM) [5], which is an extension of the Bose-Hubbard model. The 
BFHM can also be considered as a generalisation of the fermionic Hubbard model. It 
is known for the case of the Bose-Hubbard model that competition between two terms, 
one is connected with the on-site energy U and another describes the nearest-neighbour 
hopping with the tunnelling parameter t, defines the state of the system (when the 
kinetics energy dominates the ground state of the system is superfluid, in the opposite 
case the ground state is a Mott insulator) P, [7j. For the case of the BFHM phase 
diagrams are more complicated because due to the presence of fermions the effective 
interaction between bosons is generated. 

The Bose-Fermi mixtures in optical lattices have been studied using a variety of 
methods [El El HQl HH H21 [HI HH [15l HSl [17]. In [8], it was demonstrated that a two- 
dimensional mixture of bosons and fermions develops a supersolid phase (this phase 
is characterized by the simultaneous presence of a density wave and phase order in 
condensate) . The case of small fermion hopping was investigated in [9j in the framework 
of composite fermion approach and composite fermions were formed by a fermion and 
one or several bosons (bosonic holes) for attractive (repulsive) Bose-Fermi-interactions. 
In [TUl [TT] , inhomogeneous (due the the presence of the trapping potential) mixtures 
of bosons and fermions were studied. Enhancement of the superfluidity due to the 
presence of fermions was predicted in [12]. The existence of the supersolid phase was 
confirmed in [15] using quantum Monte-Carlo simulations. A mixture of the mean field 
approximation for a bosonic part and the dynamical mean field theory for a fermionic 
part of the Hamiltonian was applied in [16] and the presence of a supersolid phase 
at weak Bose-Fermi interaction was established. The case of ID Bose-Fermi-Hubbard 
model (BFHM) in the limit of large fermion hopping was investigated in [17] (the case 
of half filling was considered only and they did not observe the supersolid state). 

It should be noted that the Bose-Fermi-Hubbard-type model can also be applied 
for the description of intercalation of ions in crystals (for example, lithium intercalation 
in TiC-2 crystals). Theoretical investigation of such process in most cases were restricted 
to the numerical ab-initio and density-functional calculations [TH [191 [20] . ft was shown 
that Li is almost fully ionized once intercalated and reconstruction of electron spectrum 
at intercalation takes place. Thus, ion-electron interaction can play a significant role 
in such systems. Another interesting feature of such crystals is a displacement of the 
chemical potential at intercalation into the conduction band. As a result, these crystals 
have metallic conductivity. At intercalation of lithium in TiC-2, phase separation into Li- 
poor and Li- rich phases occurrs and this two-phase behaviour leads to a constant value 
of the electrochemical potential [2TJ, [22] (this fact is used when constructing batteries). 
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In our previous works [231 [21] we have formulated the pseudospin-electron model of 
intercalation. We have revealed that the effective interaction between bose-atoms (ions) 
can change its character depending on fermionic band filling, that leads to the charge- 
ordered phase or phase separation into the uniform phases with different concentrations 
of bosons and fermions. The ion-electron interaction was also considered in [25] at the 
investigation of thermodynamics of the spin-1 model of intercalation (the model was 
similar to the known Blume- Emery-Griffiths model), but the electron as well as ion 
transfer was not taken into account. Models of pseudospin-electron model type are 
widely used in physics of the strongly-correlated electron systems. Application of such 
models to high-temperature superconductors allows one to describe thermodynamics of 
anharmonic oxygen ion subsystem and explain the appearance of inhomogeneous states 
and the bistability phenomena (see [26]). Models of a lattice gas are also used at the 
description of ionic conductors and at the calculation of their conductivity starting from 
works of Mahan [27] and others [281129] . 

In this work we consider the hard-core limit (infinite on-site boson-boson 
interaction) of the BFHM at finite temperature (most previous investigations considered 
the zero-temperature case). Our paper is organized as follows. In section [2] we present 
the description of the model and give a self-consistent scheme for calculation of the 
density-density correlator (susceptibility) in the random phase approximation (RPA). In 
section[3]we present phase diagrams for different values of the model parameters. Special 
attention is paid to the influence of the temperature change on the phase transitions. 
We present our conclusions in section HI 



2. MODEL AND METHOD 



We consider the BFHM in the hard core limit. Using the pseudospin formalism, the 
Hamiltonian of the model is written in the following form 

= - E ^s+sr - £ Ujcfcj + gSfra 

ij ij i 

i i 

The pseudospin variable Sf takes two values (Sf = 1/2 when boson is present in a 
site i and Sf = —1/2 in the opposite case), while cf and c$ are fermionic creation 
and annihilation operators, respectively. The first and the second terms in equation 
(PQ) are responsible for the nearest neighbour boson and fermion hopping, respectively; 
g-term accounts for the boson-fermion interaction energy. To control the number of 
bosons and fermions we introduce the bosonic and fermionic chemical potentials h and 
[A, respectively. 

The unperturbed Hamiltonian Hq in the mean feald approximation (MFA) is 
obtained using the following simplification: 

gmSf -> g( ni )Sf + gn^Sf) - g^^Sf) (2) 

nstsj n{st)sj + nsf{sj) - n(st)(sj). (3) 



Phase diagrams of the Bose-Fermi-Hubbard model at finite temperature 4 

The Hamiltonain becomes 

H = H + H int , (4) 
H o = -J2 l ^ c i + E (9(S>i + 9St(n) - g(S z )(n)) 



ij 



- J^ihS? + vm) - ]T (2fyS?(S») - a,(^) 2 ) , (5) 

i ij 

i 

[(5? - - + 5?fiJ] • (6) 

u 

It is worth noting that application of the MFA to the strongly correlated systems in 
the limit of a weak on-site correlation (when there is no correlational splitting of the 
fermionic band) allows one to satisfactorily describe their properties. 

To diagonalize the Hamiltonian Hq we pass to ^-representation and perform the 
unitary transformation in the pseudospin subspace: 

S* = o* cos d + af sine, (7) 
Sf = a f cos Q-al sin9, (8) 
. „ 2Q{S X ) „ h-gn 

sin6 = X — L , cos6 = — f—, (9) 

A A 

A = v / (g(n)-h^ + (2n(S-)) 2 , Q = n q=0 , (10) 

H = - E^ fc + ^ C k C k - 22 ^ 
k i 

-Ng(S z )(n) + NQ(S x ) 2 , (11) 

where N is the number of lattice sites. 

To calculate the density-density correlator &ij{r) = (T T S* (r)Sj (0)) , we perform an 
expansion in powers of H int 



{T T Si(T)S*(0) o-((3)) 



(T T S?{T)Sm) = N r 'Y ' / (12) 



exp (-0H) = exp (-/3H )a((5), (13) 
a(/3) = T T exp - / H int {r)dr , (14) 



/ H mt (r)dr , 

(T t 5/(t)5/(0)> = <r T 5?(r)^(0)> - 

^ J dn(T T S? (r)^(0)F int (r 1 )) + .., (15) 

the averaging (...)o is performed over the distribution with Hq, where T T is the imaginary 
time ordering operator and /3 = 1/T is the inverse temperature. 

To calculate the average values of the T r -products of the pseudospin and fermion 
operators, we utilize the diagrammatic technique based on Wick's theorem for the 
spin operators [30] (besides the usual procedure for the Fermi operators). After 
elimination in this way of the nondiagonal operators we perform the semi-invariant 
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expansion in order to calculate the mean values of the remaining products of the a z 
operators. At the summation of diagrams we restrict ourselves to the diagrams having a 
structure of multi-loop chains in the spirit of the random phase approximation (see |31j). 
The junctions between bosonic (pseudospin) Green's functions and semi-invariants are 
realised by bosonic hopping Q q and the fermionic loop H q (u). It is useful to introduce 
unperturbed bosonic and fermionic Green's functions (T r cr i + (r)cr~(0))o = —2(a z )Ki m (r) 
and (T r Cfc(r)c+(0)) , respectively, and semi-invariant (T T af (t)o"^(0)) = (o~ z ) 2 + Mi m . 

Let us consider the Green's function G^f(r) = — |(T T cr"(r)(T^(0)) (with a, (3 = 
+,—,z). Typical RPA diargams for this Green's function G^(w) in the frequency 
representation are shown in figure [TJ We used the notations for the unperturbed bosonic 

K(eo) K ( ro ) ^ 




Figure 1. Typical RPA diagrams for the Green's function G^(u>). Solid and dashed 
lines with arrows denote the unperturbed bosonic and fermionic Green's functions, 
respectively. Wavy lines indicate the energy dispersion for the bosons Q q , circles and 
ovals denote the average value (er 2 ) and semi- invariants, respectively. 



Green's function 



K{u n ) = -. r, (16) 

lU) n — A 

the fermionic loop 

tt / \ 1 \- nftfc-g) - n(t k ) 

semi-invariant M(co n ) = /3<L ni o(| — (o~ z ) 2 ) and average value of the pseudospin variable 
(a z ) = |tanh(/3A/2). 

The equation for this Green's function G q ^(uj n ) is 

Gf(u n ) = G$ q (u n )A<*? + G^{u n )T%{ Un )G*(u n ), (18) 

where S^(w n ) = n^(w n ) + Qf^ (with a, (3 = +, — , z) and u n is the bosonic Matsubara 
frequency. These matrix equations (1181) form three independent sets of equations of the 
third order and can be separately solved. For the case of the Green's functions G q ~(u n ), 
G~~(ui n ), G q ~(u n ) the matrices n^(w n ), and the unperturbed Green's functions 
G$) q (un) are 

U q + (cu n ) = IL+-{u n ) = n+ + (u; n ) (19) 



Phase diagrams of the Bose- Fermi- Hubbard model at finite temperature 6 

= n— (w n ) = 3 2 n 9 k)^, (20) 

Il q z (co n ) = n+ 2 (w n ) = g 2 U q {u n ) sm9 cos 9, (21) 

nj"K) = nf K) = -/n g K)^^, (22) 

n**(w n ) = -<7 2 n g (o;„) cos 2 0, ft** = 2ft q sin 2 0, (23) 

ft q + = = -0,(1 + cos 2 0), (24) 

ft++ = ft~~ = -ft 9 (cos 2 9 - 1), (25) 

ft-* = ft+* = 2ft q sin0cos0, (26) 

ft*- = ft*+ = -ft q sin cos 0, (27) 

G+7(u; n ) = i^K)^), Gp)"K) = K(-u n )(a z ), (28) 

(4K) = MK), A+- = l, A"- = 0, A*" = 0. (29) 



Similar matrix equations can be written for the Green's functions G q + (u n ), G q + (u n ), 
G q + {uj n ) and G q z (u n ), G~ z (u n ), G z q z {u) n ) with the corresponding matrices 11^ (u n ) and 
Q q P (we do not present here these matrices). As a result, we can solve these three sets 
of equations of the third order and after some tedious algebra we derive the expression 
for the density-density correlator 

^•(r) = (T T a z (r)a z (0))cos 2 9+(T T anr)a](0))sm 2 9 

+(T t o-*(t)o-*(0)) sin # cos # + (T t <t?(t)<7?(0)) sin 9 cos G, (30) 

sin 2 9(a z ) + AM(u n ) cos 2 - 2ft 9 M(w n ) 

©q(^n) - " 

x(X-2(a z )Q q ), (31) 
D = -(iw n ) 2 + (A - 2(a z )ft q )[A - 2{a z ) cos 2 flft q 

+ (cx z )# 2 sin 2 #nq(w n ) -2M(u n )VL q \ sin 

+M(w n )( ? 2 An (? (a; n ) cos 2 9 - 2(a z )Q q M{u n )U q {u n )g 2 }. (32) 

If we use the equation of motion method developed for the two-time Green's function 
{(S z (t)\S z (t')}} = -\9{t - t')([S z (t),S z (t')\) and decoupling in the spirit of Tyablikov 
approximation [erf, H] ~ — 2i(cr 2 ) ^ erf flu + iAof we can obtain the expression for the 
correlator ((S z \S z )) q)U which is similar to the equation (l3Tj) but differs from it due 
to the absence of the terms proportional to 5 w »,o- These terms have appeared in the 
diagrammatic technique due to the presence of the semi-invariants and are important 
when we investigate the static limit u — > 0. The equation of motion method does 
not allow us to take into account these terms and because of this we should use the 
diagrammatic technique. It should be noted that such an peculiarity was also pointed 
out in [32] at the investigation of the Bose-Hubbard model in the hard-core case. 

3. PHASE DIAGRAMS 



Lines of the instability with respect to the transition into the phase with charge ordering 
can be obtained using the condition of divergence of the static density-density correlator 
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<3 q (u> = 0). We consider two cases: i) the transition from a normal (NR) nonsuperfluid 
uniform to nonsuperfluid charge- density- wave (CDW) phase ii) the transition from a 
superfluid phase to superfluid phase with long-range ordering (a supersolid phase). The 
equations for averages (n), (S z ), (S x ) are obtained in the mean field approximation. Let 
us introduce two sublattices: (n ia ) = n a , (S z a ) = (S^), a = 1,2 is a sublattice index, 
i is an elementary cell index. Using the Hamiltonian H , we can obtain the following 
equations for averages [2] 



N ^ 



1 + COs(20) ( ' hhapi 



+ 



E 

k 



(s z a ) 
(s x a ) 



1 — COs(20) f ^k£_ 

2 

h — gn c 



-i 



2X a 



tanh 



2X r 



tanh 




with 



A 



(Si) + (S 



+ 




sin 20 



(Sf)-(Sf) 



+ tl 



\a=y/(gn a -h)* + (2n{S%}y, a ±p. 
The grand canonical potential can be written as 



(33) 
(34) 
(35) 

(36) 
(37) 

(38) 



ti-\ k 



-T In 



4 cosh 



0A, 



cosh 



-g{n x (Sl) +n 2 (S*)) + 2Cl(S*)(SZ). (39) 

The doubling of the unit cell leads to the splitting in the fermionic spectrum with the gap 
g\(S$-(S*)\. The differences 5n = (th)-^), SS Z = (Sf)-(SI), 5S X = (S X )-(S%) play 
the role of the order parameter for the modulated phase ((S x ) ^ in the superfluid phase 
and SS X 7^ in the supersolid phase). Coming from the set of equations (|33|) . (1341) and 
( 13"5"j) . we can write the equations for 8n, 5S Z , 5S X and separate the contributions of the 
first order. As a result, we obtain the condition of the appearance of nonzero solutions for 
5n, 5S Z and 5S X . It can be shown that this condition coincides with the condition when 
the static density- density correlator <3 q=7T {uj = 0) diverges. Therefore our scheme for 
calculation of the density- density correlator in the RPA and the corresponding averages 
(n), (S z ) and (S x ) in the MFA is a self-consistent scheme. 
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At the numerical calculations of the density- density correlator we consider a three- 
dimensional case with a lattice constant a = 1 and in our calculations we choose a half 
width of the fermionic band W to be our energy scale (— W < tk < W). First we 
investigate the uniform phase ((ni) = (n 2 )). From the set of equations (|33|) . (134)) and 
(1351) it follows that the solutions of these equations with (S x ) 7^ can be realised when 
Q > IT . Therefore at finite temperature we can consider the transition from the normal 
uniform nonsuperfluid phase (at low temperatures this is a Mott insulating phase) to the 
CDW phase for small values of the bosonic hopping parameter (fi < 2T). In figure [2(a), 
we plot lines of the instability with respect to the transition into the charge ordered phase 
at the fixed fermionic chemical potential (the case ji = corresponds to the half-filling 
case) for the case of the nonsuperfluid phase (the bosonic concentration n B = S z + 1/2). 
As seen in figure Efa), the highest temperature of the instability is realised for the case of 
the chess-board phase with the wave vector q = (tt, tt, tt). In figure EJ^b) the lines of the 
instability for the case \x = 0.3 (when the system goes away from the half-filling case) 
are plotted. From figure [21(b) we observe that the incommensurate charge ordered phase 
with the wave vector q w (2, 2, 2) has the highest temperature of the instability and the 
system undergoes the transition to the incommensurate modulated phase. It should be 
noted that the condition of the divergence of the static density- density correlator allows 
one to investigate the phase transitions of the second order only. 



q=(Wi,7i) q=(2,2,2) 




Figure 2. The lines of instability of the nonsuperfluid phase with respect to the 
transition into the charge ordered phase for W = 1, g = —OA, O = 0, fi = (a) and 
fi = 0.3 (b). 

Now let us consider the transition to the supersolid phase. In figure [3] lines of 
instability of the superfluid phase with respect to the transition into the supersolid 
phase for the half-filling case \i = 0, ft = 0.15 are depicted. As shown in figure [31 the 
transition to the supersolid phase with modulation wave vector q = (tt, tt, tt) is realised. 
We revealed that when the system goes away from the half-filling case and /1 7^ the 
supersolid phase with the modulation wave vector q = (tt, tt, tt) also has the highest 
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temperature of the transition. 

It should be emphasised that appearance of the CDW and supersolid phases is 
connected with the presence of the effective interaction between bosons which is formed 
due to the boson-fermion correlation. This interaction depends on the filling of the 
fermionic band. 



0.05 



0.02 





q=(7r,jt,7t) 




q=(2.7,2.7,2.7) - 











0.3 0.4 0.5 0.6 0.7 

riB 



Figure 3. The lines of instability of the superfiuid phase with respect to the transition 
into the supersolid phase for W = 1, g = —0.4, fl = 0.15 and /x = 0. 

Now we want to investigate the case of the chess-board phase in more detail. We use 
the equations for averages f l33j) . ( 155]) and the expression for the grand canonical 

potential ( |39|) to find thermodynamically stable states (in this part of our numerical 
calculations we use the semielliptic density of states, pit) = -^-^y/W 2 — t 2 ). The phase 
transition lines and particle concentrations as functions of the bosonic chemical potential 
are shown in figure H] and figure [5j The phase transition from the normal uniform 
nonsuperfluid to chess-board phase can be of the second or first order, see figure Hk. 
The existence of the phase transition of the first order leads to phase separation in the 
regime of the fixed concentrations into the NR and CDW phases. As shown in figure [5j 
similar picture is obtained for the transition from the superfiuid to the supersolid phase 
and the transition from the superfiuid to the supersolid phase can be of the first or 
second order depending on the system parameters. 

In figure El we show the phase diagrams in the plane (h — Q) at low temperature. As 
temperature increases, the regions of the existence of the CDW phase and the supersolid 
phase are possible for smaller parameter space and the first order phase transitions from 
the normal uniform nonsuperfluid (superfiuid) into the CDW (SS) phases transforms 
into the second one. It should be noted that similar diagrams at T = were obtained 
in [17], but they did not reveal the possibility of the transition to the supersolid phase. 
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Figure 4. (a) The lines of the phase transitions of the second (solid lines) and hrst 
(dashed lines) order for W = 1, g = —0.4, [i = and £1 = 0. (b) The dependence of the 
bosonic and fermionic concentrations on the bosonic chemical potential at T = 0.04. 




Figure 5. (a) The lines of the phase transitions of the second (solid lines) and first 
(dashed lines) order for W = 1, g = —0.4, fj, = and £1 = 0.15. (b) The dependence of 
the bosonic and fermionic concentrations on the bosonic chemical potential at T = 0.02. 



4. CONCLUSIONS 



The phase transitions in the Bose-Fermi-Hubbard model at finite temperature has been 
considered in this work. We studied the hard-core limit and used pseudospin formalism. 
The thermodynamics of the model was investigated in the case of the weak boson- 
fermion interaction. The analytical expression for the density-density correlator has 
been obtained in the framework of the self-consistent scheme of the random phase 
approximation. The effective boson-boson interaction is formed due to the boson 
interaction with fermions, this effective interaction depends on the filling of the fermionic 
band. It is revealed that at small values of the bosonic tunnelling amplitude the system 
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Figure 6. Phase diagram in the (h — ft) plane for W = 1, g = —0.4, /i = and 
T = 0.01. Solid (dashed) lines denote the phase transitions of the second (first) order. 

undergoes the phase transition from the uniform nonsuperfluid phase to the chess-board 
phase (the case of half filling of the fermion band) or to the incommensurate phase 
(when the system goes away from the half filling) at the lowering of the temperature. At 
increase of the bosonic hopping parameter the phase transition from the superfluid phase 
to the supersolid phase with a doubly modulated lattice period takes place (it should be 
noted that the presence of the supersolid phase at the weak boson-fermion interaction 
and zero-temperature was also established in [16] in the framework of a generalized 
dynamical mean field theory). The transition from the uniform to modulated phase can 
be of the first or second order, depending on the model parameters and temperature. 
The presence of the first order phase transition means that in the regime of the fixed 
fermionic concentrations the phase separation into the uniform and modulated phases 
is possible. 
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